% get three circle shape from an image
im = imread('threeCircles.png');

% inside == 1, outside == 0
dout = bwdist( im ); din = -bwdist( 1-im );
% put them together to one distance map
dall = zeros( size( im ) );
dall( im==1 ) = din( im==1 );
dall( im==0 ) = dout( im==0 );

nrIter = 200;
dt = 0.1;
a = -1;

phi = dall; z = zeros( size( im ) );
phixP = z; phixM = z; phiyP = z; phiyM = z; phiX = z; phiY = z;

for i=0:nrIter
if ( mod(i,20)==0 )
figure, contour( phi,[0 0], 'k', 'linewidth', 2 ); axis equal; box on;
set( gca, 'xticklabel', '' ); set( gca, 'yticklabel', '');
end
% compute the derivatives (forward P and backward M)
phixP(1:end-1,:) = phi(2:end,:)-phi(1:end-1,:);
phixM(2:end,:) = phi(2:end,:)-phi(1:end-1,:);
phiyP(:,1:end-1) = phi(:,2:end)-phi(:,1:end-1);
phiyM(:,2:end) = phi(:,2:end)-phi(:,1:end-1);
% get the indices and derivatives for x direction
indxPP = find( a*phixP>0 & a*phixM>0 );
indxMM = find( a*phixP<0 & a*phixM<0 );
indxZ = find( a*phixM<=0 & a*phixP>=0 );
indxS = find( a*phixM>=0 & a*phixP<=0 );
indxSGP = find( abs(phixP(indxS))>=abs(phixM(indxS)) );
indxSGM = find( abs(phixM(indxS))>=abs(phixP(indxS)) );
phix( indxPP ) = phixM( indxPP );
phix( indxMM ) = phixP( indxMM );
phix( indxZ ) = 0;
phix( indxS( indxSGP ) ) = phixP( indxS( indxSGP ) );
phix( indxS( indxSGM ) ) = phixM( indxS( indxSGM ) );
% get the indices and derivatives for y direction
indxPP = find( a*phiyP>0 & a*phiyM>0 );
indxMM = find( a*phiyP<0 & a*phiyM<0 );
indxZ = find( a*phiyM<=0 & a*phiyP>=0 );
indxS = find( a*phiyM>=0 & a*phiyP<=0 );
indxSGP = find( abs(phiyP(indxS))>=abs(phiyM(indxS)) );
indxSGM = find( abs(phiyM(indxS))>=abs(phiyP(indxS)) );
phiy( indxPP ) = phiyM( indxPP );
phiy( indxMM ) = phiyP( indxMM );
phiy( indxZ ) = 0;
phiy( indxS( indxSGP ) ) = phiyP( indxS( indxSGP ) );
phiy( indxS( indxSGM ) ) = phiyM( indxS( indxSGM ) );
% now do the update
phi = phi - a*dt*( (phix.^2 + phiy.^2).^0.5 );
end